Biomolecular interactions determine how transcription factors recognize their DNA binding sites, how proteins interact with each other, and consequently how a biological system functions. Since many biological molecules bear considerable electric charge, electrostatic interactions are among the most important when studying biomolecular interactions. However, electrostatic interactions in biological systems are difficult to calculate accurately in practice. Aside from the significants charges carried by biomolecules such as DNA and proteins, the solvent itself namely, water produces considerable electrostatic effects. Furthermore, hydrogen bonds, known to be involved in helix formation in both DNA and proteins, are essentially electrostatic in origin. Indeed, it seems that electrostatic effects often drive the physical-chemical processes in biological systems and, thereby, determine biological function. Therefore, any attempt to perform molecular dynamics (MD) simulations of biological systems will require an adequate description of these electrostatic forces. In the past, we developed a controllable scheme for computing the electrostatic energy. In principle, one may take the numerical derivative to obtain the electrostatic force. In reality, the numerical derivative computations are expensive and may become even more expensive when a high degree of accuracy is demanded. In order to have a practical system for simulation of biomolecular systems, it becomes necessary to calculate the electrostatic forces by function evaluation of analytical formulae, which is one of the primary efforts of this fiscal year. In the process of doing this, we have also implemented a more stable method for evaluating mathematical functions in a rotated coordinate system. We have also developed a new approach that allows for exact calculation of energetics even if the dielectric constant changes in space in a continuous fashion. This scenario is applicable to the the context, for example, of water layer around biomolecules. To the best of our knowledge, our method is the only one that uses energetic functional in a throughout consistent manner. We have writen up this result and the paper was published in Physical Review E.